* ==============================================================================
* ==============================================================================
* Title: Evaulating the effects of IMF conditionality: An extension of quantitative approaches and an empirical application to public education spending
* Authors: Thomas Stubbs, Alexander Kentikelenis, Bernhard Reinsberg, Lawrence King
* Version: Final (16 October 2018)
* ==============================================================================
* ==============================================================================
* STATA code for Appendix A of 'Evaluating the effects of IMF conditionality'
/*
  Cross-combination of two parameters
  1. Severity of bias  
  -- instrument weakness
  -- cross correlation
  2. Omitted variables?
  -- no unobserved confounder
  -- unobserved confounders in one stage (predictor X4)
  
  Performance of different models
  a. CFA-IV/MLE
  d. IV-IV/MLE
  b. CFA-/MLE
  e. -IV/2SLS
  i. IV-/2SLS
  f. OLS
*/

  set matsize 800
  set seed 12345
  cmp setup
  clear

* *********************************************************************
* Setup 1
* *********************************************************************

* Create matrix to store results 
  mat A=J(1,14,0)
  mat D=A
  mat B=J(1,10,0)
  mat E=B
  mat I=B
  mat F=J(1,8,0)
  
  mat cr1=J(500,2,0)
  mat F1=J(500,6,0)
  
  forvalues sims = 1/500 {
  clear
  qui set obs 1000

  qui g id=.
  forvalues i=1/50{
  local a=1+20*(`i'-1)
  local b=20*`i'
  qui replace id=`i' in `a'/`b'
  } 

  forvalues i=1/50{
  qui g id`i'=(id==`i')
  }
  qui drop id

* Specify parameters
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  
  
* Specify the correlation and generate variables
  local y13=0.2
  local y23=0.2
  local y12=0.2
  matrix C = (1,`y12',`y13' \ `y12',1,`y23' \ `y13',`y23',1)
  drawnorm u1 u2 u3, means(0,0,0) sds(1,1,1) corr(C)

  gen X1=rnormal(0,1)
  gen X2=rnormal(0,1)
  gen X3=rnormal(0,1)
  gen e1=rnormal(0,1)
  gen e2=rnormal(0,1)
  
  local eta1=0.3
  local eta2=0.3
  
  gen Z1=`eta1'*X1+(1-`eta1')*e1
  qui reg X1 Z1 
  mat Coef=e(b)
  mat cr1[`sims',1]=Coef[1,1]^2
  
  gen Z2=`eta2'*X2+(1-`eta2')*e2
  qui reg X2 Z2
  mat Coef=e(b)
  mat cr1[`sims',2]=Coef[1,1]^2
  

* Generate responses
  gen y1star=b10+b11*X1 +u1  -id5 -id3 -id2 +id17 +id36 +id48
  egen y1starmedian=median(y1star)
  gen y1=(y1star>y1starmedian)
  drop y1starmedian
  gen y2=b20+b21*X2+ u2 -2.4*id1 -2.3*id2 -2.2*id3 -2.1*id4 -2.0*id5 -1.9*id6 -1.8*id7 -1.7*id8 -1.6*id9 -1.5*id10 -1.4*id11 -1.3*id12 -1.2*id13 -1.1*id14 -1*id15 -.9*id16 -.8*id17 -.7*id18 -.6*id19 -.5*id20 -.4*id21 -.3*id22 -.2*id23 -.1*id24 +.1*id26 +.2*id27 +.3*id28 +.4*id29 +.5*id30 +.6*id31 +.7*id32 +.8*id33 +.9*id34 +id35 +1.1*id36 +1.2*id37 +1.3*id38 +1.4*id39 +1.5*id40 +1.6*id41 +1.7*id42 +1.8*id43 +1.9*id44 +2*id45 +2.1*id46 +2.2*id47 +2.3*id48 +2.4*id49 +2.5*id50
  replace y2=0 if y1==0
  gen y3=b30+b31*y1 +b32*y2 +b33*X3 +u3 +1.1*id1 +1.8*id2 +2.4*id3 +1.6*id4 -.4*id5 +1.7*id6 -1.1*id7 -.3*id8 +.5*id9 +1.5*id10 -2.3*id11 -1.9*id12 +.3*id13 -2.4*id14 +2.3*id15 -.9*id16 -1.2*id17 +.8*id18 +1.3*id19 +.4*id20 +id21 -1.5*id22 +.2*id23 -.6*id24 -.8*id25 +1.2*id26 -1.3*id28 -2*id29 -.1*id30 +.9*id31 +2.2*id32 -2.1*id33 -1.4*id34 +.6*id35 -id36 +1.7*id37 +.1*id38 -1.8*id39 +1.4*id40 +1.9*id41 -.2*id42 +.7*id43 +2*id44 -.5*id45 +2.1*id46 -.7*id47 -2.2*id48 -1.6*id49

  
* Run the models saving beta estimates and SEs and store results

  ** CFA-IV/MLE
  qui cmp (y1=Z1 X3) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_probit $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F1[`sims',1]=r(chi2)
  mat A=(A \ betas, vces)
  
  
  ** IV-IV/MLE
  qui cmp (y1=Z1 X3 id*) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_cont $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F1[`sims',2]=r(chi2)
  mat D=(D \ betas, vces)
  
  
  ** CFA
  qui treatreg y3 y1 y2 X3 id*, treat(y1=Z1 X3) twostep hazard(imr)
  qui reg y3 y1 y2 X3 imr id*
  drop imr
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[imr]
  mat betas[1,5]=_b[_cons]
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[imr]^2
  mat vces[1,5]=_se[_cons]^2
  
  mat B=(B \ betas, vces)
  
  
  ** 2SLS: conditions
  qui cmp (y3=y1 y2 X3 id*) (y2=Z2 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z2
  mat F1[`sims',4]=r(chi2)
    
  mat E=(E \ betas, vces)
  
  
  ** 2SLS: programs
  qui cmp (y3=y1 y2 X3 id*) (y1=Z1 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z1
  mat F1[`sims',5]=r(chi2)
    
  mat I=(I \ betas, vces)
  
  
  ** OLS
  qui reg y3 y1 y2 X3 id*
  mat betas=J(1,4,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[_cons]
  mat vces=J(1,4,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[_cons]^2
  mat F=(F \ betas, vces)
  
  di `sims'
}

  mat li F1
  mat li cr1

* Results for A
  mat A = A[2...,.]
  clear 
  svmat A, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A1: read off new line here (only the 'mean' column is relevant)
  su MBias* RMSE* opt*
  
  
* Results for D
  mat D= D[2...,.]
  clear 
  svmat D, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*

  ** Table A1: read off new line here (only the 'mean' column is relevant)
  su MBias* RMSE* opt*

  
* Results for B
  mat B=B[2...,.]
  clear 
  svmat B, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A1: read off new line here (only the 'mean' column is relevant)
  su MBias* RMSE* opt*

  
* Results for I
  mat I=I[2...,.]
  clear 
  svmat I, names(col)
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A1: read off new line here (only the 'mean' column is relevant)
  su MBias* RMSE* opt*
  
  
* Results for E
  mat E=E[2...,.]
  clear 
  svmat E, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A1: read off new line here (only the 'mean' column is relevant)
  su MBias* RMSE* opt*
  
  
* Results for OLS
  mat F=F[2...,.]
  clear 
  svmat F, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c5)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c6)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*

  ** Table A1: read off new line here (only the 'mean' column is relevant)
  su MBias* RMSE* opt*  
	
	
* *********************************************************************
* Setup 2
* *********************************************************************
	
* Create matrix to store results 
  mat A=J(1,14,0)
  mat D=A
  mat B=J(1,10,0)
  mat E=B
  mat I=B
  mat F=J(1,8,0)
  
  mat cr2=J(500,2,0)
  mat F2=J(500,6,0)
  
  forvalues sims = 1/500 {
  clear
  qui set obs 1000

  qui g id=.
  forvalues i=1/50{
  local a=1+20*(`i'-1)
  local b=20*`i'
  qui replace id=`i' in `a'/`b'
  } 

  forvalues i=1/50{
  qui g id`i'=(id==`i')
  }
  qui drop id

* Specify parameters
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  
  
* Specify the correlation and generate variables
  local y13=0.2
  local y23=0.2
  local y12=0.2
  matrix C = (1,`y12',`y13' \ `y12',1,`y23' \ `y13',`y23',1)
  drawnorm u1 u2 u3, means(0,0,0) sds(1,1,1) corr(C)

  gen X1=rnormal(0,1)
  gen X2=rnormal(0,1)
  gen X3=rnormal(0,1)
  gen e1=rnormal(0,1)
  gen e2=rnormal(0,1)
  
  local eta1=0.3
  local eta2=0.08
  
  gen Z1=`eta1'*X1+(1-`eta1')*e1
  qui reg X1 Z1 
  mat Coef=e(b)
  mat cr2[`sims',1]=Coef[1,1]^2
  
  gen Z2=`eta2'*X2+(1-`eta2')*e2
  qui reg X2 Z2
  mat Coef=e(b)
  mat cr2[`sims',2]=Coef[1,1]^2
  

* Generate responses
  gen y1star=b10+b11*X1 +u1  -id5 -id3 -id2 +id17 +id36 +id48
  egen y1starmedian=median(y1star)
  gen y1=(y1star>y1starmedian)
  drop y1starmedian
  gen y2=b20+b21*X2+ u2 -2.4*id1 -2.3*id2 -2.2*id3 -2.1*id4 -2.0*id5 -1.9*id6 -1.8*id7 -1.7*id8 -1.6*id9 -1.5*id10 -1.4*id11 -1.3*id12 -1.2*id13 -1.1*id14 -1*id15 -.9*id16 -.8*id17 -.7*id18 -.6*id19 -.5*id20 -.4*id21 -.3*id22 -.2*id23 -.1*id24 +.1*id26 +.2*id27 +.3*id28 +.4*id29 +.5*id30 +.6*id31 +.7*id32 +.8*id33 +.9*id34 +id35 +1.1*id36 +1.2*id37 +1.3*id38 +1.4*id39 +1.5*id40 +1.6*id41 +1.7*id42 +1.8*id43 +1.9*id44 +2*id45 +2.1*id46 +2.2*id47 +2.3*id48 +2.4*id49 +2.5*id50
  replace y2=0 if y1==0
  gen y3=b30+b31*y1 +b32*y2 +b33*X3 +u3 +1.1*id1 +1.8*id2 +2.4*id3 +1.6*id4 -.4*id5 +1.7*id6 -1.1*id7 -.3*id8 +.5*id9 +1.5*id10 -2.3*id11 -1.9*id12 +.3*id13 -2.4*id14 +2.3*id15 -.9*id16 -1.2*id17 +.8*id18 +1.3*id19 +.4*id20 +id21 -1.5*id22 +.2*id23 -.6*id24 -.8*id25 +1.2*id26 -1.3*id28 -2*id29 -.1*id30 +.9*id31 +2.2*id32 -2.1*id33 -1.4*id34 +.6*id35 -id36 +1.7*id37 +.1*id38 -1.8*id39 +1.4*id40 +1.9*id41 -.2*id42 +.7*id43 +2*id44 -.5*id45 +2.1*id46 -.7*id47 -2.2*id48 -1.6*id49
	
  ** CFA-IV/MLE
  qui cmp (y1=Z1 X3) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_probit $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F2[`sims',1]=r(chi2)
  mat A=(A \ betas, vces)
  
  
  ** IV-IV/MLE
  qui cmp (y1=Z1 X3 id*) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_cont $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F2[`sims',2]=r(chi2)
  mat D=(D \ betas, vces)
  
  
  ** CFA
  qui treatreg y3 y1 y2 X3 id*, treat(y1=Z1 X3) twostep hazard(imr)
  qui reg y3 y1 y2 X3 imr id*
  drop imr
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[imr]
  mat betas[1,5]=_b[_cons]
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[imr]^2
  mat vces[1,5]=_se[_cons]^2
  
  mat B=(B \ betas, vces)
  
  
  ** 2SLS: conditions
  qui cmp (y3=y1 y2 X3 id*) (y2=Z2 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z2
  mat F2[`sims',4]=r(chi2)
    
  mat E=(E \ betas, vces)
  
  
  ** 2SLS: programs
  qui cmp (y3=y1 y2 X3 id*) (y1=Z1 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z1
  mat F2[`sims',5]=r(chi2)
    
  mat I=(I \ betas, vces)
  
  
  ** OLS
  qui reg y3 y1 y2 X3 id*
  mat betas=J(1,4,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[_cons]
  mat vces=J(1,4,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[_cons]^2
  mat F=(F \ betas, vces)
  
  di `sims'
}

  mat li F2
  mat li cr2

* Results for A
  mat A = A[2...,.]
  clear 
  svmat A, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*

  ** Table A2--new line  
  su MBias* RMSE* opt*
  
  
* Results for D
  mat D= D[2...,.]
  clear 
  svmat D, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*

  ** Table A2--new line  
  su MBias* RMSE* opt*

  
* Results for B
  mat B=B[2...,.]
  clear 
  svmat B, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A2--new line
  su MBias* RMSE* opt*
  
* Results for I
  mat I=I[2...,.]
  clear 
  svmat I, names(col)
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A2--new line
  su MBias* RMSE* opt*
  
  
* Results for E
  mat E=E[2...,.]
  clear 
  svmat E, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A2--new line
  su MBias* RMSE* opt*
    
  
* Results for OLS
  mat F=F[2...,.]
  clear 
  svmat F, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c5)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c6)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A2--new line
  su MBias* RMSE* opt*
    
	
		
	
* *********************************************************************
* Setup 3
* *********************************************************************

* Create matrix to store results 
  mat A=J(1,14,0)
  mat D=A
  mat B=J(1,10,0)
  mat E=B
  mat I=B
  mat F=J(1,8,0)
  
  mat cr3=J(500,2,0)
  mat F3=J(500,6,0)
  
  forvalues sims = 1/500 {
  clear
  qui set obs 1000

  qui g id=.
  forvalues i=1/50{
  local a=1+20*(`i'-1)
  local b=20*`i'
  qui replace id=`i' in `a'/`b'
  } 

  forvalues i=1/50{
  qui g id`i'=(id==`i')
  }
  qui drop id

* Specify parameters
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  
  
* Specify the correlation and generate variables
  local y13=0.2
  local y23=0.2
  local y12=0.2
  matrix C = (1,`y12',`y13' \ `y12',1,`y23' \ `y13',`y23',1)
  drawnorm u1 u2 u3, means(0,0,0) sds(1,1,1) corr(C)

  gen X1=rnormal(0,1)
  gen X2=rnormal(0,1)
  gen X3=rnormal(0,1)
  gen e1=rnormal(0,1)
  gen e2=rnormal(0,1)
  
  local eta1=0.08
  local eta2=0.08
  
  gen Z1=`eta1'*X1+(1-`eta1')*e1
  qui reg X1 Z1 
  mat Coef=e(b)
  mat cr3[`sims',1]=Coef[1,1]^2
  
  gen Z2=`eta2'*X2+(1-`eta2')*e2
  qui reg X2 Z2
  mat Coef=e(b)
  mat cr3[`sims',2]=Coef[1,1]^2
  

* Generate responses
  gen y1star=b10+b11*X1 +u1  -id5 -id3 -id2 +id17 +id36 +id48
  egen y1starmedian=median(y1star)
  gen y1=(y1star>y1starmedian)
  drop y1starmedian
  gen y2=b20+b21*X2+ u2 -2.4*id1 -2.3*id2 -2.2*id3 -2.1*id4 -2.0*id5 -1.9*id6 -1.8*id7 -1.7*id8 -1.6*id9 -1.5*id10 -1.4*id11 -1.3*id12 -1.2*id13 -1.1*id14 -1*id15 -.9*id16 -.8*id17 -.7*id18 -.6*id19 -.5*id20 -.4*id21 -.3*id22 -.2*id23 -.1*id24 +.1*id26 +.2*id27 +.3*id28 +.4*id29 +.5*id30 +.6*id31 +.7*id32 +.8*id33 +.9*id34 +id35 +1.1*id36 +1.2*id37 +1.3*id38 +1.4*id39 +1.5*id40 +1.6*id41 +1.7*id42 +1.8*id43 +1.9*id44 +2*id45 +2.1*id46 +2.2*id47 +2.3*id48 +2.4*id49 +2.5*id50
  replace y2=0 if y1==0
  gen y3=b30+b31*y1 +b32*y2 +b33*X3 +u3 +1.1*id1 +1.8*id2 +2.4*id3 +1.6*id4 -.4*id5 +1.7*id6 -1.1*id7 -.3*id8 +.5*id9 +1.5*id10 -2.3*id11 -1.9*id12 +.3*id13 -2.4*id14 +2.3*id15 -.9*id16 -1.2*id17 +.8*id18 +1.3*id19 +.4*id20 +id21 -1.5*id22 +.2*id23 -.6*id24 -.8*id25 +1.2*id26 -1.3*id28 -2*id29 -.1*id30 +.9*id31 +2.2*id32 -2.1*id33 -1.4*id34 +.6*id35 -id36 +1.7*id37 +.1*id38 -1.8*id39 +1.4*id40 +1.9*id41 -.2*id42 +.7*id43 +2*id44 -.5*id45 +2.1*id46 -.7*id47 -2.2*id48 -1.6*id49

  ** CFA-IV/MLE
  qui cmp (y1=Z1 X3) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_probit $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F3[`sims',1]=r(chi2)
  mat A=(A \ betas, vces)
  
  
  ** IV-IV/MLE
  qui cmp (y1=Z1 X3 id*) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_cont $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F3[`sims',2]=r(chi2)
  mat D=(D \ betas, vces)
  
  
  ** CFA
  qui treatreg y3 y1 y2 X3 id*, treat(y1=Z1 X3) twostep hazard(imr)
  qui reg y3 y1 y2 X3 imr id*
  drop imr
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[imr]
  mat betas[1,5]=_b[_cons]
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[imr]^2
  mat vces[1,5]=_se[_cons]^2
  
  mat B=(B \ betas, vces)
  
  
  ** 2SLS: conditions
  qui cmp (y3=y1 y2 X3 id*) (y2=Z2 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z2
  mat F3[`sims',4]=r(chi2)
    
  mat E=(E \ betas, vces)
  
  
  ** 2SLS: programs
  qui cmp (y3=y1 y2 X3 id*) (y1=Z1 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z1
  mat F3[`sims',5]=r(chi2)
    
  mat I=(I \ betas, vces)
  
  
  ** OLS
  qui reg y3 y1 y2 X3 id*
  mat betas=J(1,4,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[_cons]
  mat vces=J(1,4,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[_cons]^2
  mat F=(F \ betas, vces)
  
  di `sims'
}

  mat li F3
  mat li cr3

* Results for A
  mat A = A[2...,.]
  clear 
  svmat A, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A3--new line
  su MBias* RMSE* opt*
  
  
* Results for D
  mat D= D[2...,.]
  clear 
  svmat D, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*

  ** Table A3--new line  
  su MBias* RMSE* opt*

  
* Results for B
  mat B=B[2...,.]
  clear 
  svmat B, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A3--new line
  su MBias* RMSE* opt*
  
* Results for I
  mat I=I[2...,.]
  clear 
  svmat I, names(col)
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A3--new line
  su MBias* RMSE* opt*
  
  
* Results for E
  mat E=E[2...,.]
  clear 
  svmat E, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*

  ** Table A3--new line
  su MBias* RMSE* opt*
    
  
* Results for OLS
  mat F=F[2...,.]
  clear 
  svmat F, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c5)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c6)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*

  ** Table A3--new line  
  su MBias* RMSE* opt*
    	

	
* *********************************************************************
* Setup 4
* *********************************************************************

* Create matrix to store results 
  mat A=J(1,14,0)
  mat D=A
  mat B=J(1,10,0)
  mat E=B
  mat I=B
  mat F=J(1,8,0)
  
  mat cr4=J(500,2,0)
  mat F4=J(500,6,0)
  
  forvalues sims = 1/500 {
  clear
  qui set obs 1000

  qui g id=.
  forvalues i=1/50{
  local a=1+20*(`i'-1)
  local b=20*`i'
  qui replace id=`i' in `a'/`b'
  } 

  forvalues i=1/50{
  qui g id`i'=(id==`i')
  }
  qui drop id

* Specify parameters
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  
  
* Specify the correlation and generate variables
  local y13=0.5
  local y23=0.6
  local y12=0.7
  matrix C = (1,`y12',`y13' \ `y12',1,`y23' \ `y13',`y23',1)
  drawnorm u1 u2 u3, means(0,0,0) sds(1,1,1) corr(C)

  gen X1=rnormal(0,1)
  gen X2=rnormal(0,1)
  gen X3=rnormal(0,1)
  gen e1=rnormal(0,1)
  gen e2=rnormal(0,1)
  
  local eta1=0.3
  local eta2=0.3
  
  gen Z1=`eta1'*X1+(1-`eta1')*e1
  qui reg X1 Z1 
  mat Coef=e(b)
  mat cr4[`sims',1]=Coef[1,1]^2
  
  gen Z2=`eta2'*X2+(1-`eta2')*e2
  qui reg X2 Z2
  mat Coef=e(b)
  mat cr4[`sims',2]=Coef[1,1]^2
  

* Generate responses
  gen y1star=b10+b11*X1 +u1  -id5 -id3 -id2 +id17 +id36 +id48
  egen y1starmedian=median(y1star)
  gen y1=(y1star>y1starmedian)
  drop y1starmedian
  gen y2=b20+b21*X2+ u2 -2.4*id1 -2.3*id2 -2.2*id3 -2.1*id4 -2.0*id5 -1.9*id6 -1.8*id7 -1.7*id8 -1.6*id9 -1.5*id10 -1.4*id11 -1.3*id12 -1.2*id13 -1.1*id14 -1*id15 -.9*id16 -.8*id17 -.7*id18 -.6*id19 -.5*id20 -.4*id21 -.3*id22 -.2*id23 -.1*id24 +.1*id26 +.2*id27 +.3*id28 +.4*id29 +.5*id30 +.6*id31 +.7*id32 +.8*id33 +.9*id34 +id35 +1.1*id36 +1.2*id37 +1.3*id38 +1.4*id39 +1.5*id40 +1.6*id41 +1.7*id42 +1.8*id43 +1.9*id44 +2*id45 +2.1*id46 +2.2*id47 +2.3*id48 +2.4*id49 +2.5*id50
  replace y2=0 if y1==0
  gen y3=b30+b31*y1 +b32*y2 +b33*X3 +u3 +1.1*id1 +1.8*id2 +2.4*id3 +1.6*id4 -.4*id5 +1.7*id6 -1.1*id7 -.3*id8 +.5*id9 +1.5*id10 -2.3*id11 -1.9*id12 +.3*id13 -2.4*id14 +2.3*id15 -.9*id16 -1.2*id17 +.8*id18 +1.3*id19 +.4*id20 +id21 -1.5*id22 +.2*id23 -.6*id24 -.8*id25 +1.2*id26 -1.3*id28 -2*id29 -.1*id30 +.9*id31 +2.2*id32 -2.1*id33 -1.4*id34 +.6*id35 -id36 +1.7*id37 +.1*id38 -1.8*id39 +1.4*id40 +1.9*id41 -.2*id42 +.7*id43 +2*id44 -.5*id45 +2.1*id46 -.7*id47 -2.2*id48 -1.6*id49

  ** CFA-IV/MLE
  qui cmp (y1=Z1 X3) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_probit $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F4[`sims',1]=r(chi2)
  mat A=(A \ betas, vces)
  
  
  ** IV-IV/MLE
  qui cmp (y1=Z1 X3 id*) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_cont $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F4[`sims',2]=r(chi2)
  mat D=(D \ betas, vces)
  
  
  ** CFA
  qui treatreg y3 y1 y2 X3 id*, treat(y1=Z1 X3) twostep hazard(imr)
  qui reg y3 y1 y2 X3 imr id*
  drop imr
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[imr]
  mat betas[1,5]=_b[_cons]
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[imr]^2
  mat vces[1,5]=_se[_cons]^2
  
  mat B=(B \ betas, vces)
  
  
  ** 2SLS: conditions
  qui cmp (y3=y1 y2 X3 id*) (y2=Z2 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z2
  mat F4[`sims',4]=r(chi2)
    
  mat E=(E \ betas, vces)
  
  
  ** 2SLS: programs
  qui cmp (y3=y1 y2 X3 id*) (y1=Z1 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z1
  mat F4[`sims',5]=r(chi2)
    
  mat I=(I \ betas, vces)
  
  
  ** OLS
  qui reg y3 y1 y2 X3 id*
  mat betas=J(1,4,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[_cons]
  mat vces=J(1,4,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[_cons]^2
  mat F=(F \ betas, vces)
  
  di `sims'
}

  mat li F4
  mat li cr4

* Results for A
  mat A = A[2...,.]
  clear 
  svmat A, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*

  ** Table A4--new line  
  su MBias* RMSE* opt*
  
  
* Results for D
  mat D= D[2...,.]
  clear 
  svmat D, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A4--new line
  su MBias* RMSE* opt*

  
* Results for B
  mat B=B[2...,.]
  clear 
  svmat B, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A4--new line
  su MBias* RMSE* opt*
  
* Results for I
  mat I=I[2...,.]
  clear 
  svmat I, names(col)
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A4--new line
  su MBias* RMSE* opt*
  
  
* Results for E
  mat E=E[2...,.]
  clear 
  svmat E, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A4--new line
  su MBias* RMSE* opt*
    
  
* Results for OLS
  mat F=F[2...,.]
  clear 
  svmat F, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c5)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c6)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A4--new line
  su MBias* RMSE* opt*
    

	
* *********************************************************************
* Setup 5 -- misspecified model due to Cor(A,C)>0
* *********************************************************************

* Create matrix to store results 
  mat A=J(1,14,0)
  mat D=A
  mat B=J(1,10,0)
  mat E=B
  mat I=B
  mat F=J(1,8,0)
  
  mat cr5=J(500,2,0)
  mat F5=J(500,6,0)
  
  forvalues sims = 1/500 {
  clear
  qui set obs 1000

  qui g id=.
  forvalues i=1/50{
  local a=1+20*(`i'-1)
  local b=20*`i'
  qui replace id=`i' in `a'/`b'
  } 

  forvalues i=1/50{
  qui g id`i'=(id==`i')
  }
  qui drop id

* Specify parameters
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  
  
* Specify the correlation and generate variables
  local y13=0.2
  local y23=0.2
  local y12=0.2
  matrix C = (1,`y12',`y13' \ `y12',1,`y23' \ `y13',`y23',1)
  drawnorm u1 u2 u3, means(0,0,0) sds(1,1,1) corr(C)

  gen X1=rnormal(0,1)
  gen X2=rnormal(0,1)
  gen X3=rnormal(0,1)
  gen e1=rnormal(0,1)
  gen e2=rnormal(0,1)
  gen X4=rnormal(0,1)
  
  local eta1=0.3
  local eta2=0.3
  
  gen Z1=`eta1'*X1+(1-`eta1')*e1
  qui reg X1 Z1 
  mat Coef=e(b)
  mat cr5[`sims',1]=Coef[1,1]^2
  
  gen Z2=`eta2'*X2+(1-`eta2')*e2
  qui reg X2 Z2
  mat Coef=e(b)
  mat cr5[`sims',2]=Coef[1,1]^2
  

* Generate responses
  gen y1star=2*X4+ b10+b11*X1 +u1  -id5 -id3 -id2 +id17 +id36 +id48
  egen y1starmedian=median(y1star)
  gen y1=(y1star>y1starmedian)
  drop y1starmedian
  gen y2=b20+b21*X2+ u2 -2.4*id1 -2.3*id2 -2.2*id3 -2.1*id4 -2.0*id5 -1.9*id6 -1.8*id7 -1.7*id8 -1.6*id9 -1.5*id10 -1.4*id11 -1.3*id12 -1.2*id13 -1.1*id14 -1*id15 -.9*id16 -.8*id17 -.7*id18 -.6*id19 -.5*id20 -.4*id21 -.3*id22 -.2*id23 -.1*id24 +.1*id26 +.2*id27 +.3*id28 +.4*id29 +.5*id30 +.6*id31 +.7*id32 +.8*id33 +.9*id34 +id35 +1.1*id36 +1.2*id37 +1.3*id38 +1.4*id39 +1.5*id40 +1.6*id41 +1.7*id42 +1.8*id43 +1.9*id44 +2*id45 +2.1*id46 +2.2*id47 +2.3*id48 +2.4*id49 +2.5*id50
  replace y2=0 if y1==0
  gen y3=2*X4+ b30+b31*y1 +b32*y2 +b33*X3 +u3 +1.1*id1 +1.8*id2 +2.4*id3 +1.6*id4 -.4*id5 +1.7*id6 -1.1*id7 -.3*id8 +.5*id9 +1.5*id10 -2.3*id11 -1.9*id12 +.3*id13 -2.4*id14 +2.3*id15 -.9*id16 -1.2*id17 +.8*id18 +1.3*id19 +.4*id20 +id21 -1.5*id22 +.2*id23 -.6*id24 -.8*id25 +1.2*id26 -1.3*id28 -2*id29 -.1*id30 +.9*id31 +2.2*id32 -2.1*id33 -1.4*id34 +.6*id35 -id36 +1.7*id37 +.1*id38 -1.8*id39 +1.4*id40 +1.9*id41 -.2*id42 +.7*id43 +2*id44 -.5*id45 +2.1*id46 -.7*id47 -2.2*id48 -1.6*id49


  ** CFA-IV/MLE
  qui cmp (y1=Z1 X3) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_probit $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F5[`sims',1]=r(chi2)
  mat A=(A \ betas, vces)
  
  
  ** IV-IV/MLE
  qui cmp (y1=Z1 X3 id*) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_cont $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F5[`sims',2]=r(chi2)
  mat D=(D \ betas, vces)
  
  
  ** CFA
  qui treatreg y3 y1 y2 X3 id*, treat(y1=Z1 X3) twostep hazard(imr)
  qui reg y3 y1 y2 X3 imr id*
  drop imr
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[imr]
  mat betas[1,5]=_b[_cons]
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[imr]^2
  mat vces[1,5]=_se[_cons]^2
  
  mat B=(B \ betas, vces)
  
  
  ** 2SLS: conditions
  qui cmp (y3=y1 y2 X3 id*) (y2=Z2 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z2
  mat F5[`sims',4]=r(chi2)
    
  mat E=(E \ betas, vces)
  
  
  ** 2SLS: programs
  qui cmp (y3=y1 y2 X3 id*) (y1=Z1 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z1
  mat F5[`sims',5]=r(chi2)
    
  mat I=(I \ betas, vces)
  
  
  ** OLS
  qui reg y3 y1 y2 X3 id*
  mat betas=J(1,4,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[_cons]
  mat vces=J(1,4,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[_cons]^2
  mat F=(F \ betas, vces)
  
  di `sims'
}

  mat li F5
  mat li cr5

* Results for A
  mat A = A[2...,.]
  clear 
  svmat A, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A5--new line
  su MBias* RMSE* opt*
  
  
* Results for D
  mat D= D[2...,.]
  clear 
  svmat D, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*

  ** Table A5--new line
  su MBias* RMSE* opt*

  
* Results for B
  mat B=B[2...,.]
  clear 
  svmat B, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A5--new line
  su MBias* RMSE* opt*
  
* Results for I
  mat I=I[2...,.]
  clear 
  svmat I, names(col)
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A5--new line
  su MBias* RMSE* opt*
  
  
* Results for E
  mat E=E[2...,.]
  clear 
  svmat E, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A5--new line
  su MBias* RMSE* opt*
    
  
* Results for OLS
  mat F=F[2...,.]
  clear 
  svmat F, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c5)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c6)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A5--new line
  su MBias* RMSE* opt*
    



* *********************************************************************
* Setup 6 -- misspecified model excess BC-correlation 
* *********************************************************************

* Create matrix to store results 
  mat A=J(1,14,0)
  mat D=A
  mat B=J(1,10,0)
  mat E=B
  mat I=B
  mat F=J(1,8,0)
  
  mat cr6=J(500,2,0)
  mat F6=J(500,6,0)
  
  forvalues sims = 1/500 {
  clear
  qui set obs 1000

  qui g id=.
  forvalues i=1/50{
  local a=1+20*(`i'-1)
  local b=20*`i'
  qui replace id=`i' in `a'/`b'
  } 

  forvalues i=1/50{
  qui g id`i'=(id==`i')
  }
  qui drop id

* Specify parameters
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  
  
* Specify the correlation and generate variables
  local y13=0.2
  local y23=0.2
  local y12=0.2
  matrix C = (1,`y12',`y13' \ `y12',1,`y23' \ `y13',`y23',1)
  drawnorm u1 u2 u3, means(0,0,0) sds(1,1,1) corr(C)

  gen X1=rnormal(0,1)
  gen X2=rnormal(0,1)
  gen X3=rnormal(0,1)
  gen e1=rnormal(0,1)
  gen e2=rnormal(0,1)
  gen X4=rnormal(0,1)
  
  local eta1=0.3
  local eta2=0.3
  
  gen Z1=`eta1'*X1+(1-`eta1')*e1
  qui reg X1 Z1 
  mat Coef=e(b)
  mat cr6[`sims',1]=Coef[1,1]^2
  
  gen Z2=`eta2'*X2+(1-`eta2')*e2
  qui reg X2 Z2
  mat Coef=e(b)
  mat cr6[`sims',2]=Coef[1,1]^2
  

* Generate responses
  gen y1star= b10+b11*X1 +u1  -id5 -id3 -id2 +id17 +id36 +id48
  egen y1starmedian=median(y1star)
  gen y1=(y1star>y1starmedian)
  drop y1starmedian
  gen y2=2*X4+ b20+b21*X2+ u2 -2.4*id1 -2.3*id2 -2.2*id3 -2.1*id4 -2.0*id5 -1.9*id6 -1.8*id7 -1.7*id8 -1.6*id9 -1.5*id10 -1.4*id11 -1.3*id12 -1.2*id13 -1.1*id14 -1*id15 -.9*id16 -.8*id17 -.7*id18 -.6*id19 -.5*id20 -.4*id21 -.3*id22 -.2*id23 -.1*id24 +.1*id26 +.2*id27 +.3*id28 +.4*id29 +.5*id30 +.6*id31 +.7*id32 +.8*id33 +.9*id34 +id35 +1.1*id36 +1.2*id37 +1.3*id38 +1.4*id39 +1.5*id40 +1.6*id41 +1.7*id42 +1.8*id43 +1.9*id44 +2*id45 +2.1*id46 +2.2*id47 +2.3*id48 +2.4*id49 +2.5*id50
  replace y2=0 if y1==0
  gen y3=2*X4+ b30+b31*y1 +b32*y2 +b33*X3 +u3 +1.1*id1 +1.8*id2 +2.4*id3 +1.6*id4 -.4*id5 +1.7*id6 -1.1*id7 -.3*id8 +.5*id9 +1.5*id10 -2.3*id11 -1.9*id12 +.3*id13 -2.4*id14 +2.3*id15 -.9*id16 -1.2*id17 +.8*id18 +1.3*id19 +.4*id20 +id21 -1.5*id22 +.2*id23 -.6*id24 -.8*id25 +1.2*id26 -1.3*id28 -2*id29 -.1*id30 +.9*id31 +2.2*id32 -2.1*id33 -1.4*id34 +.6*id35 -id36 +1.7*id37 +.1*id38 -1.8*id39 +1.4*id40 +1.9*id41 -.2*id42 +.7*id43 +2*id44 -.5*id45 +2.1*id46 -.7*id47 -2.2*id48 -1.6*id49

  ** CFA-IV/MLE
  qui cmp (y1=Z1 X3) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_probit $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F6[`sims',1]=r(chi2)
  mat A=(A \ betas, vces)
  
  
  ** IV-IV/MLE
  qui cmp (y1=Z1 X3 id*) (y2=Z2 X3 id*) (y3=y1 y2 X3 id*), indicators($cmp_cont $cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,7,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_2:_cons]
  mat betas[1,4]=_b[lnsig_3:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  mat betas[1,6]=_b[atanhrho_13:_cons]  
  mat betas[1,7]=_b[atanhrho_23:_cons]
  
  mat vces=J(1,7,0) 
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_2:_cons]^2
  mat vces[1,4]=_se[lnsig_3:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2
  mat vces[1,6]=_se[atanhrho_13:_cons]^2  
  mat vces[1,7]=_se[atanhrho_23:_cons]^2
  
  qui test (Z1 Z2)
  mat F6[`sims',2]=r(chi2)
  mat D=(D \ betas, vces)
  
  
  ** CFA
  qui treatreg y3 y1 y2 X3 id*, treat(y1=Z1 X3) twostep hazard(imr)
  qui reg y3 y1 y2 X3 imr id*
  drop imr
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[imr]
  mat betas[1,5]=_b[_cons]
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[imr]^2
  mat vces[1,5]=_se[_cons]^2
  
  mat B=(B \ betas, vces)
  
  
  ** 2SLS: conditions
  qui cmp (y3=y1 y2 X3 id*) (y2=Z2 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z2
  mat F6[`sims',4]=r(chi2)
    
  mat E=(E \ betas, vces)
  
  
  ** 2SLS: programs
  qui cmp (y3=y1 y2 X3 id*) (y1=Z1 X3 id*), indicators($cmp_cont $cmp_cont) iterate(50)
  mat betas=J(1,5,0)
  mat betas[1,1]=_b[y3:y1]
  mat betas[1,2]=_b[y3:y2]
  mat betas[1,3]=_b[lnsig_1:_cons]
  mat betas[1,4]=_b[lnsig_2:_cons]
  mat betas[1,5]=_b[atanhrho_12:_cons] 
  
  mat vces=J(1,5,0)
  mat vces[1,1]=_se[y3:y1]^2
  mat vces[1,2]=_se[y3:y2]^2
  mat vces[1,3]=_se[lnsig_1:_cons]^2
  mat vces[1,4]=_se[lnsig_2:_cons]^2
  mat vces[1,5]=_se[atanhrho_12:_cons]^2 
  qui test Z1
  mat F6[`sims',5]=r(chi2)
    
  mat I=(I \ betas, vces)
  
  
  ** OLS
  qui reg y3 y1 y2 X3 id*
  mat betas=J(1,4,0)
  mat betas[1,1]=_b[y1]
  mat betas[1,2]=_b[y2]
  mat betas[1,3]=_b[X3]
  mat betas[1,4]=_b[_cons]
  mat vces=J(1,4,0)
  mat vces[1,1]=_se[y1]^2
  mat vces[1,2]=_se[y2]^2
  mat vces[1,3]=_se[X3]^2
  mat vces[1,4]=_se[_cons]^2
  mat F=(F \ betas, vces)
  
  di `sims'
}

  mat li F6
  mat li cr6

* Results for A
  mat A = A[2...,.]
  clear 
  svmat A, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A6--new line
  su MBias* RMSE* opt*
  
  
* Results for D
  mat D= D[2...,.]
  clear 
  svmat D, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c8)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c9)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A6--new line
  su MBias* RMSE* opt*

  
* Results for B
  mat B=B[2...,.]
  clear 
  svmat B, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A6--new line
  su MBias* RMSE* opt*
  
* Results for I
  mat I=I[2...,.]
  clear 
  svmat I, names(col)
  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A6--new line
  su MBias* RMSE* opt*
  
  
* Results for E
  mat E=E[2...,.]
  clear 
  svmat E, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c6)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c7)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A6--new line
  su MBias* RMSE* opt*
    
  
* Results for OLS
  mat F=F[2...,.]
  clear 
  svmat F, names(col)

  g b10=0.5
  g b11=1
  
  g b20=-0.5
  g b21=2 
  
  g b30=1
  g b31=1.5
  g b32=-0.5
  g b33=-3
  note: closing state `c(seed)'

  ** Bias
  gen Bias1 = 1/(c1/b31)
  gen Bias2 = 1/(c2/b32)
  egen MBias1=mean(Bias1)
  egen MBias2=mean(Bias2)
  drop Bias*

  ** RMSE
  gen SqEr1=(c1-b31)^2
  gen SqEr2=(c2-b32)^2
  egen MSE1 = mean(SqEr1)
  egen MSE2 = mean(SqEr2)
  gen RMSE1 = sqrt(MSE1)
  gen RMSE2 = sqrt(MSE2)
  drop SqEr* MSE*
  
  ** Optimism
  egen mean1=mean(c1)
  gen SqEr1=(c1-mean1)^2
  egen TotSqEr1=total(SqEr1)
  gen OptTop1=sqrt(TotSqEr1)
  egen sumvar1=total(c5)
  gen OptBot1=sqrt(sumvar1)
  gen opt1=OptTop1/OptBot1
  
  egen mean2=mean(c2)
  gen SqEr2=(c2-mean2)^2
  egen TotSqEr2=total(SqEr2)
  gen OptTop2=sqrt(TotSqEr2)
  egen sumvar2=total(c6)
  gen OptBot2=sqrt(sumvar2)
  gen opt2=OptTop2/OptBot2
  drop mean* SqEr* Tot* sumvar* Opt*
  
  ** Table A6--new line
  su MBias* RMSE* opt*
    